Ultracold atoms at unitarity within quantum Monte Carlo 
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Variational and diffusion quantum Monte Carlo (VMC and DMC) calculations of the properties of the zero- 
temperature fermionic gas at unitarity are reported. The ratio of the energy of the interacting to the non- 
interacting gas for a system of 128 particles is calculated to be 0.4517(3) in VMC and 0.4339(1) in the more 
accurate DMC method. The spherically-averaged pair-correlation functions, momentum densities, and one- 
body density matrices are very similar in VMC and DMC, but the two-body density matrices and condensate 
fractions show some differences. Our best estimate of the condensate fraction of 0.51 is a little smaller than 
values from other quantum Monte Carlo calculations. 
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I. INTRODUCTION 

Superfluid pairing in ultracold trapped atoms has been the 
subject of much experimental and theoretical work.i^iMl 
The range of the inter-atomic interaction in a dilute atomic 
Fermi gas is much smaller than the average distance between 
the atoms, and only the i-wave scattering length a is rele- 
vant. The only relevant dimensionless coupling parameter is 
then l/{akp), where kp is the Fermi wave vector of the gas. 
The scattering length can be altered by applying a magnetic 
field and, by using Fano-Feshbach resonanceS j^i^i^i^ 1 / (akp) 
may be varied from large negative to large positive values. 
When the interaction is weakly attractive and kp is sufficiently 
small, 1 / {akp) — > — oo, and the gas is in the Bardeen-Cooper- 
Schrieffer (BCS) superfluid regime. When I /{akp) 
the molecules are tightly bound and the system forms a Bose- 
Einstein condensate (BEC). The behavior in the intermedi- 
ate regime changes smoothly with I /{akp) and at unitarity, 
where |l/(a^F)| ^ 0, a smooth crossover between BCS-Uke 
and BEC-like behavior occurs. 

At unitarity the scattering length becomes larger than the 
inter-particle distance and the only energy scale is kp (we con- 
sider a particle mass of unity). The ground state energy Eq 
can therefore be written as 



^0 -^—kl, 



(1) 



where the factor of 3/10 is chosen so that ^ is the fraction of 
the energy of the non-interacting Fermi gas at the same den- 
sity. A number of experimental and theoretical determinations 
of the universal parameter ^ have been reported. In each case 
the parameter was found to be smaller than unity, showing that 
the interactions are attractive at unitarity. 

In this paper we report calculations of the energy, pair cor- 
relation functions (PCFs), momentum density and the one- 
and two-body density matrices of the Fermi gas at unitar- 
ity. We use the zero-temperature variational and diffusion 



quantum Monte Carlo methodsi^ (VMC and DMC), as have 
been used in previous studies i''^''^'''*''^''^i'^''^ Our study dif- 
fers from earlier ones mainly in the construction of the trial 
wave functions, the larger system size used, and in studying 
the dependence of the energy on the particle density. Other 
quantum Monte Carlo (QMC) methods have been used to 
study ultracold atomic systems at finite temperatures.'^-^" 

The rest of this paper is set out as follows. In Sec.|ll]we dis- 
cuss the Hamiltonian used to model the atomic Fermi gas. An 
introduction to the VMC and DMC methods is given in Sec. 
nil A! and specific points pertaining to Fermi atomic gases are 
described in Sec. lIIIBl Our trial wave functions are described 
in Sec. IIIICI and important parameters of the DMC calcula- 
tions are discussed in Sec. IIIIDI Our results are reported and 
discussed in Sec. lIVI and the conclusions of our study are sum- 
marized in Sec.lv] 



II. THE HAMILTONIAN 



The Hamiltonian takes the form 



(2) 



where v(ry ) is the interaction potential. We use face-centered 
cubic (fee) simulation cells subject to periodic boundary con- 
ditions. We wish to study the system with a delta-function 
potential, but this is difficult to sample using Monte Carlo 
methods. We have instead used the Poschl-Teller interaction 
which, on resonance, is given by 



v{nj) 



cosh (/ir,j 



(3) 
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where r,y is the distance between particles / and j, and 2//i 
is the effective width of the potential well. Since the inter- 
particle interaction is very short-ranged, particles of the same 
spin are kept apart by the antisymmetry of the wave func- 
tion, and the interaction between them is negligible for the 
well widths used here and would be precisely zero for the 
delta-function potential. We therefore set the interaction be- 
tween particles of the same spin to zero, as has been done 
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in previous calculations. The Poschl-Teller interaction has 
been used in previous QMC calculations/^''^ and we pre- 
fer it to the square-well which has also been used in QMC 
calculations/^''^ because its smoothness aids Monte Carlo 
sampling. In all of our QMC calculations reported here we 
have used jJ. — 12. 

The particle density is kp/{37i^), but we report densities in 
terms of the rj parameter, which is the radius of a sphere con- 
taining one atom on average, andrj — {97t/4-y^^/kp. For most 
of our calculations we have used a density parameter = 1, 
so that jXTs ^ 1, as required for dilute conditions, although in 
Sec. IIV Dl we report some investigations of the effect of in- 
creasing r,.. 



where £l(R) = ^-j-'/Z^t is the local energy, and p = py in 
VMC and p^p^ in DMC. 

DMC expectation values of operators which do not com- 
mute with the Hamiltonian depend on the entire trial wave 
function, not just its nodal surface. To reduce this bias, at the 
expense of increasing the noise, one can use the extrapolation 
approximation,— 

(^) ^ 2 (A)^^^ - (A) VMC + ^ [(^T - O)'] , (7) 

where (4)omc and (^)vmc DMC and VMC ex- 

pectation values of operator A, respectively. The quantity 
( )dmc ^ ( )vMC gi^^^ ^ measure of the accuracy of the trial 
wave function. 



III. QMC METHODS 

A. VMC and DMC methods 

In VMC the energy is calculated as the expectation value of 
the Hamiltonian with an approximate many-body trial wave 
function. In the more accurate DMC method the ground state 
energy is obtained by evolving the wave function in imagi- 
nary time so that it decays towards the ground state. Projector 
methods such as DMC suffer from a fermion sign problem, 
which is evaded by making the fixed-node approximation, 
and importance sampling is introduced to reduce the statisti- 
cal noise. The importance-sampled fixed-node fermion DMC 
algorithm was first used by Ceperley and Alder to study the 
electron gas.^' 

The key quantity in VMC and DMC calculations is the 
trial wave function, which controls both the accuracy that 
can be obtained and the statistical efficiency of the compu- 
tation. In VMC the accuracy of the energy estimate depends 
on the whole of the trial wave function, while in DMC it de- 
pends only on the form of its nodal surface, as the DMC algo- 
rithm (in principle) gives the lowest energy compatible with 
the fixed nodal surface. In practice the DMC energy estimate 
also shows some dependence on the timestep used and a very 
weak dependence on the population of particle configurations. 
Improving the trial wave function tends to reduce these biases 
and improve the statistical efficiency of the calculations. 

The VMC algorithm generates particle configurations dis- 
tributed according to 



/'v(R) = 



(R) 



J^'l{R')dR' ' 



(4) 



where is a real trial wave function and R is the 3A^- 
dimensional vector of the coordinates of the particles. DMC 
generates configurations distributed according to 



Pd{R) = 



^t(R)'I'(R) 

7*iVR0*(RVR' ' 



(5) 



where ^(R) is the DMC wave function. The total energy in 
both VMC and DMC is calculated from 



p(R)£l(R)^R, 



(6) 



B. VMC and DMC calculations for Fermi atomic gases 

The construction of accurate trial wave functions for Fermi 
atomic gases is not straightforward. The variation of the wave 
function must be described very accurately at small inter- 
particle separations where the interaction potential varies very 
rapidly. The binding energy of an isolated molecule is vanish- 
ingly small at resonance, but the exact value of ^ for the gas is 
certainly smaller than the BCS mean-field value of 0.59, and 
therefore the interactions between molecules are very impor- 
tant at unitarity. The exact wave function for an isolated pair 
of opposite spin fermionic atoms at resonance decays as the 
inverse distance between the particles, and we must describe 
the deviations from this behavior in the gas phase. We con- 
clude that it is necessary to provide a good description of both 
the long- and short-ranged behavior of the pairing function to 
obtain accurate results for the system at unitarity. 

The simulations are performed with a finite number of par- 
ticles, and we wish to obtain results which accurately reflect 
those that would be obtained with an infinite number. A great 
deal of experience has been gained in performing QMC calcu- 
lations for electron gases, and the finite size effects are greatly 
reduced by averaging the energies obtained at different wave 
vectors^^ ("twist averaging") and extrapolating the averaged 
energies to the infinite system limit. The QMC studies of ul- 
tracold atoms reported so far have not employed very large 
numbers of particles and have not used twist averaging, and 
it is not clear whether the results are converged with respect 
to system size. Although we have not used twist averaging 
in the calculations reported here, we have used systems with 
128 particles, which is approximately twice the largest num- 
ber used in previous DMC calculations, '^''^''''''^''^''^''** 

Another problem is that we really want the solution for a 
delta-function potential, but for computational reasons we use 
a well of finite width. The ground state of the many-particle 
system with a delta-function potential is a molecular gas be- 
cause bound states with more than two particles cannot exist 
in this case but, with a finite well-width, clusters of particles 
can form at high densities. In practice this instability of the 
gas phase occurs only for densities where jj.rs ^ 1, and we 
work at much lower densities. Nevertheless, it is clear that 
the results of QMC calculations will depend on both the well 
width and the particle density. 
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C. Trial wave functions 

We used a singlet-pairing BCS-like wave function consist- 
ing of a determinant of identical pairing orbitals, (p (r,y), each 
of which is a function of the separation of an up- and a down- 
spin particle. The pairing orbital is represented by a sum of 
polynomial terms, 

We set (p to zero for rij greater than the cutoff length Lp. The 
third power in Eq. ^ was chosen to ensure that the local en- 
ergy is continuous at Lp. The 7,, are optimizable parameters, 
but the value of 71 is determined by the condition that (p is 
cuspless at the origin. We tested various values of A^p and 
chose A^p = 4 for the results presented here. The optimized 
value of the cutoff length of Lp = 3.5 a.u. is 3.5 times the 
average distance between particles. We also tested orbitals 
represented by linear combinations of Gaussian orbitals, lin- 
ear combinations of plane waves, and linear combinations of 
Gaussian orbitals, plane waves, and polynomials. Gaussian 
orbitals also appeared to be a useful choice, but we chose the 
polynomial basis as the orbitals and their derivatives can be 
evaluated much more rapidly. 

The determinant of pairing orbitals is multiplied by a Jas- 
trow factor of the form e'', with 

-/M = L^Lexp{iG.r,-,}+U-^ ^0„r;;, (9) 

s=l GGS \ i-J / „=o 

where A^s is the number of stars of symmetry-related recipro- 
cal lattice vectors G, and Nj is the order of the polynomial. 
The Xs and 9„ are optimizable parameters, although 9i is de- 
termined by the condition that J is cuspless at the origin, and 
we set the polynomial part of J to zero for r,,- > Lj.^^ After 
some testing we chose Ns = 4- an N] — 8, and we used an op- 
timized value of L] = 0.86 a.u. 

We also applied backflow a transformation to the determi- 
nant of pairing orbitals. The particle coordinates r, are 
replaced by collective coordinates x,(R) = r, + ^i{R), where 
^,(R) is the backflow displacement of particle /, which de- 
pends on all the particle positions. The backflow displacement 
is given by 

C,(R) = I%(n;)ro'- (10) 
We have used the form 

= (^)'Lp«'-r,-, (11) 

where Lg is a cutoff length and A^b is the order of the poly- 
nomial in the backflow expansion and the p„ are optimizable 
parameters, with pi determined by the condition that 77 is cus- 
pless at the origin. We chose Nb —5 and used an optimized 
value of Lb = 1 04 a.u. 



The wave functions were optimized within a VMC proce- 
dure by minimizing the mean absolute deviation of the lo- 
cal energies from their median value. We found this ap- 
proach to be superior to variance minimization schemes. 
The polynomial term in the Jastrow factor is much more im- 
portant than the plane-wave part. The Jastrow and backflow 
functions can, in principle, operate between both parallel and 
anti-parallel spin particles, although the correlation effects be- 
tween the parallel-spin particles beyond the exchange interac- 
tion already included in the determinant are small. We did, 
however, find a small lowering of the VMC energy when we 
allowed the plane wave parameters in the Jastrow factor to 
be non-zero for parallel-spin particles. We did not include 
parallel-spin terms for the polynomials. The wave function 
contains a total of 28 variable parameters. 

D. QMC calculations 

We used the CASINO coda^ for all of our QMC calcu- 
lations. We performed some test calculations with 32 and 64 
particles, but all of our results reported in this paper were ob- 
tained with 128 particles. We used a time step of 0.015 a.u. for 
all the DMC results presented in this paper. Test calculations 
using a timestep of 0.03 a.u. did not change the total energy 
within the statistical error bars achieved. We used a mean 
population of 3200 configurations, and tests indicated that the 
population control bias with this number of configurations is 
negligible. 

IV. RESULTS 
A. Total energy and the £, parameter 

When evaluating the £, parameter it is not immediately ob- 
vious whether to use the non-interacting energy Lni of the 
finite system studied, the infinite system, or some other value. 
Energies of non-interacting systems for various particle num- 
bers and fee and simple cubic (sc) simulation cells are shown 
in Fig. [T] Lni oscillates in an irregular manner about the 
infinite system value as the particle number A^ is increased 
and converges rather slowly with A^. Earlier DMC calcula- 
tions of ^ used sc cells and particle numbers from A^ = 14 
to 6642ii3ji6ji7 rpj^g error in L^i for the 66-particle system is 
quite small (0.5%), although it increases to nearly 5% for the 
sc cell with 80 particles. The rate of convergence for the fee 
cell is similar to that for the sc cell. 

The convergence of the interacting energy with system size 
has not been well-studied, but using more particles is expected 
to improve the results. The momentum distribution of the non- 
interacting system has a discontinuity at the Fermi momen- 
tum, while that of the interacting system is smooth, see Sec. 
llVCi so that the kinetic energy of the non-interacting system 
varies rapidly in k-space, leading to rapid fluctuations in the 
kinetic energy with system size. It therefore seems likely that 
the interacting energy will converge faster with system size 
than the non-interacting energy. 
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FIG. 1 : (Color online) Energy per particle of the non-interacting sys- 
tem at r, = 1 for different particle numbers, for the sc and fee lattices. 
The dotted line at 1.10495 a.u. shows the exact energy of the infinite 
system. 



Method 




Exp.^ 


0.32(10) 


ExpiS 


0.51(4) 


Exp,^ 


0.46(5) 


Theori^i 


0.455 


Theor,^ 


0.360(20) 


DMC12 


0.44(1) 


DMCM 


0.42(1) 


DMGi^ 


0.414(5) 


DMC'^ 


0.42(1) 


VMC This work 


0.4517(3) 


DMC This work 


0.4339(1) 



TABLE I: Values of the universal parameter £, from experiments and 
theory. The error bars given in brackets for the current work are 
purely statistical errors, and they do not account for biases from finite 
size effects or fixed-node errors or other sources. 



Our values of B,, and those from other calculations and 
some experiments, are given in Table U The DMC energy 
is bounded from above by the VMC energy and from below 
by the exact energy and, as expected, the VMC values of ^ 
are a little larger than the DMC ones. The results reported 
in Table |l] were obtained using the value of ^ni for the fi- 
nite system studied of 1.08307 a.u. per particle, compared 
with the infinite-system result of 1.10495 a.u. per particle. 
As discussed above, one might argue that it is more appro- 
priate to calculate ^ using the infinite-system value of ^ni, as 
the interacting energy is expected to converge faster with sys- 
tem size than the non-interacting energy. Using the infinite- 
system non-interacting energy gives values of ^ — 0.4428(1) 
in VMC, and t, = 0.4253(1) in DMC, which are about 2% 
lower than those reported in Table U Our DMC value of ^ is 
very similar to previous ones. ^^'^^''^''^ 




FIG. 2: Raw data for the parallel-spin PCF. 




FIG. 3: Raw data for the anti-parallel-spin PCF. 



B. Pair correlation functions 

We evaluated the spatially and rotationally averaged pair 
correlation functions (PCFs) for the parallel and anti-parallel 
spin pairs, which are shown in Figs. |2] and [3] The difference 
in our VMC and DMC results was negligible on the scales 
of Figs. 12] and [3] and consequently the effect of using the ex- 
trapolation of Eq. (|7]l is negligible. Figures |2] and [3] show the 
DMC data itself, not fits to the data. The noise in the data 
is very small, but can just be resolved in Fig.|2]at small r/r,, 
where the number of counts is small. The parallel-spin PCF 
shown in Fig. |2] shows a hole largely confined to the region 
r/rs < 2, which is essentially an exchange hole, and the PCF is 
very similar to the non-interacting (Hartree-Fock) result. This 
is consistent with the fact that we found only a very weak 
parallel-spin Jastrow factor The anti-parallel-spin PCF (Fig. 
|3]l shows a very strong enhancement for r/r, < 1 arising from 
the pairing. The behavior at small r/r, is not shown as it de- 
pends strongly on the well width. The anti-parallel PCF dips 
below unity in the region 1 < r/r^ < 2. The PCFs are similar 
to those reported in Fig. 3 of Astrakharchik et al. and Fig. 1 
of Chang and Pandharipande,J^ 
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C. Momentum density and density matrices 



The one-body density-matxix (OBDM) may be written as 



dr2 . . . dri\ 



Jp{R)dR 

and the two-body density-matrix (TBDM) as 



(12) 



p3(ri,r2;r'i, 



dvT, . . . dr^ 



!p{K)dR 



(13) 

where p{K) is the VMC or DMC probabiHty disti'ibution of 
Eqs. (IHl or (|5]l, ri and r'j are a-spin particle coordinates, r2 
and r'2 are j3 -spin particle coordinates, andA^a =Np =N/2is 
the number of particles of each spin type. 

We have evaluated the translationally and rotationally av- 
eraged density matrices, which we denote by Pa'(r) and 
(2) 

Pa {r), respectively. The momentum density n{k) is the 

Fourier transform of Pa\r), but we evaluate it directly in 
Fourier space, which is a somewhat better numerical ap- 
proach. Our data for n{k) at unitarity are broadly similar to 
those presented in Fig. 2 of Ref. 12 and Fig. 2 of Ref. 15, al- 
though in detail there are some differences. Our data show 
a monotonic decrease of n [k) with increasing momentum, in 
common with the results of Ref. ^T3, but in conflict with Ref. 
[l2i which show a peak below the Fermi momentum. Carlson 
et alP show VMC and DMC data for 14 and 38 particles. 
The VMC and DMC data are in good agreement above k-p, but 
below k-p there are substantial differences. The differences be- 
tween our VMC and DMC data are very small, and would be 
barely visible on the scale of Fig.|4] This suggests that our trial 
wave functions are superior to those of Ref.ll^. Astrakharchik 
et al}^ only report DMC data in the form of a fit, rather than 
giving the calculated values. Our momentum density is a little 
closer to the BCS form than the curve shown in Ref. 

The OBDM is shown in Fig.|5] Again, the VMC and DMC 
data are virtually indistinguishable so that extrapolation is un- 
necessary. Our calculated OBDM is very similar to that shown 
in Fig. 1 of Ref. [H 

The condensate fraction c is related to the translationally 
and rotationally averaged TBDM by 



(14) 



where Q. is the volume of the simulation cell and N/2 is the 
number of pairs of particles in the system. VMC and DMC 
data for the TBDM are shown in Fig. |6] In this case we find 
some differences between the VMC and DMC results. We ob- 
tained a condensate fraction of c = 0.57 in VMC and c = 0.54 
in DMC, see Fig. |6] so that the value obtained using the ex- 
trapolation of Eq. (|7]) is 0.5 1 . Our DMC values are lower than 
those of 0.61(2) for 38 particles and 0.57(2) for 66 particles 
reported by Astrakharchik et alJ^ 




FIG. 4: (Color online) Momentum density. The raw DMC data are 
shown (crosses) and the dotted line is a fit to this data. 




FIG. 5: Raw DMC data for the one-body density-matrix. 
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FIG. 6: (Color online) The two-body density-matrix. The VMC data 
are shown as plus signs and the DMC data as squares. The horizon- 
tal lines shows estimates of the asymptotic behavior which give the 
condensate fraction. 
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D. Varying tlie particle density 

As mentioned in Sec. [Ill we require jirs ^ 1 for dilute con- 
ditions. This can be satisfied by, for example, fixing r, and 
choosing the effective width of the potential well l/jx to be 
sufficiently small, or by fixing l/jj. and choosing to be suf- 
ficiently large. We have tried both of these approaches, but 
did not obtain a smooth variation of the energy when reduc- 
ing the range of the interaction, at least partly because the 
wave function optimization becomes more difficult. We ob- 
tained smoother results when increasing the value of while 
keeping jJ. — 12, as shown in Table HH and the optimizations 
worked well in these calculations. The value of £, slowly in- 
creases with jUrj, but this behavior could arise from the fixed- 
node error inherent in the trial wave functions. Note that the 
differences between the VMC and DMC energies increase 
with jxrs, indicating that the trial wave functions are becom- 
ing less accurate. Increasing the size of the simulation cell 
extends the range of the trial wave function and makes it more 
difficult to represent. 





S, (VMC) E, (DMC) 


12 


0.456(1) 


0.4370(4) 


14 


0.462(1) 


0.4379(4) 


16 


0.470(1) 


0.4392(6) 


18 


0.477(1) 


0.4442(5) 



TABLE II: The energy parameter, with fl = 12 and different par- 
ticle density parameters r^. 



teraction at unitarity using 128 particles, which is larger than 
in previous calculations. Our DMC result of ^ = 0.4339(1) 
is similar to previous DMC results, li-i^ but is significantly 
larger than that of ^ = 0.360(20) obtained from a recent ap- 
plication of the epsilon expansion. The VMC and DMC re- 
sults for the spherically-averaged pair-correlation functions, 
the momentum density, and the one-body density matrix were 
in good agreement, illustrating the high accuracy of our trial 
wave functions. The VMC and DMC results for the two-body 
density matrix, and the condensate fraction derived from it, 
are somewhat different, indicating that significant errors still 
arise from the approximate trial wave functions and/or finite 
size effects. We have calculated a somewhat smaller conden- 
sate fraction than in other studies using similar methods. We 
also calculated the variation of ^ with particle density for a 
fixed well width, finding a relatively small variation over the 
range of densities studied, but we were unable to draw a firm 
conclusion as to whether these represent real variations with 
well width or whether they are due to variations in the fixed- 
node errors. 

Acknowledgments 



V. CONCLUSIONS 

We have performed VMC and DMC calculations of the 
atomic Fermi gas at zero temperature with a short ranged in- 



We are grateful to Neil Drummond for fruitful discussions. 
This work was supported by the Engineering and Physical Sci- 
ences Research Council (EPSRC) of the UK. Computational 
resources were provided by the Cambridge High Performance 
Computing Service. 



1 W. Ketterle, Rev. Mod. Phys. 74, 1131 (2002). 

2 E. A. Cornell and C. E. Wieman, Rev. Mod. Phys. 74, 875 (2002). 

3 T. Kohler, K. Goral, and R S. lulienne. Rev. Mod. Phys. 78, 1311 
(2006). 

^ S. Giorgini, L. P. Pitaevskii, and S. Stringari, Rev. Mod. Phys. 80, 
1215 (2008). 

5 I. Bloch, I. Dalibard, and W. Zwerger, Rev. Mod. Phys. 80, 885 
(2008). 

^ U. Fano, Phys. Rev. 124, 1866 (1961). 

^ H. Feshbach, Ann. Phys. (N.Y) 19, 287 (1962). 

^ S. Inouye, M. R. Andrews, J. Stenger, H.-J. Miesner, D. M. 

Stamper-Kurn, andW. Ketterle, Nature 392, 151 (1998). 
^ P. Courteille, R. S. Freeland, D. I. Heinzen, F. A. van Abeelen, 

and B. I. Verhaar, Phys. Rev. Lett. 81, 69 (1998). 

We use Hartree atomic units, so that energies are in Ha and masses 

are in multiples of the electron mass. 
" W. M. C. Foulkes, L. Mitas, R. I. Needs, and G. Rajagopal, Rev. 

Mod. Phys. 73, 33 (2001). 

I. Carlson, S.-Y. Chang, V. R. Pandharipande, and K. E. Schmidt, 
Phys. Rev. Lett. 91, 050401 (2003). 



'■^ G. E. Astrakharchik, I. Boronat, I. Casulleras, and S. Giorgini, 
Phys. Rev. Lett. 93, 200404 (2004). 

S. Y. Chang, V. R. Pandharipande, I. Carlson, and K. E. Schmidt, 
Phys. Rev. A 70, 043602 (2004). 

G. E. Astrakharchik, I. Boronat, I. Casulleras, and S. Giorgini, 
Phys. Rev. Lett. 95, 230405 (2005). 
1'' S. Y. Chang and V. R. Pandharipande, Phys Rev. Lett. 95, 080402 
(2005). 

1^ I. Carlson and S. Reddy, Phys. Rev. Lett. 95, 060401 (2005). 
A. Gezerlis, S. Gandolfi, K. E. Schmidt, and I. Carlson, Phys. 
Rev. Lett. 103, 060403 (2009). 

V. K. Akkineni, D. M. Ceperley, and N. Trivedi, Phys. Rev. B 76, 
165116 (2007). 

^" E.Burovski, E. Kozik, N Prokof 'ev, B. Svistunov, and M. Troyer, 

Phys. Rev. Lett. 101, 090402 (2008). 
2' D. M. Ceperley and B. I. Alder, Phys. Rev. Lett. 45, 566 (1980). 

D. M. Ceperley and M. H. Kalos, in Monte Carlo Methods in 

Statistical Physics, 2nd ed., edited by K. Binder (Springer, Berlin, 

Germany, 1986), p. 145. 
23 C. Lin, F H. Zong, and D. M. Ceperley, Phys. Rev. E 64, 016702 



7 



(2001). 

N. D. Dmmmond, M. D. Towler, and R. J. Needs, Phys. Rev. B 
70, 235119 (2004). 

R. R Feynman, Phys. Rev. 94, 262 (1954). 

R. R Feynman and M. Cohen, Phys. Rev. 102, 1189 (1956). 

P. Lopez Rios, A. Ma, N. D. Drummond, M. D. Towler, and R. J. 

Needs, Phys. Rev. E 74, 066701 (2006). 

R R. C. Kent, R. J. Needs, and G. Rajagopal, Phys. Rev. B 59, 

12344 (1999). 

N. D. Drummond and R. J. Needs, Phys. Rev. B 72, 085124 
(2005). 

R. J. Needs, M. D. Towler, N. D. Drunmiond, and P. Lopez Rfos, 



CASINO version 2.4 User Manual, University of Cambridge, 
Cambridge (2009). 

M. Bartenstein, A. Altmeyer, S. Riedl, S. Jochim, C. Chin, J. H. 
Denschlag, and R. Grimm, Phys. Rev. Lett. 92, 120401 (2004). 
J. Kinast, A. Turlapov, J. E. Thomas, Q. J. Chen, J. Stajic, and K. 
Levin, Science 307, 1296 (2005). 

G. B. Partridge, W. Li, R. I. Kamar, Y. Liao, and R. G. Hulet, 
Science 311, 503 (2006). 

A. Perali, R Fieri, and G. C. Strinati, Phys. Rev. Lett. 93, 100404 
(2004). 

Y. Nishida, Phys. Rev. A 79, 013627 (2009). 



